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^ . Abstract 

We consider a chain of one-dimensional dipole moments connected to two thermal baths with 
^ ' different temperatures. The system is in nonequilibrium steady state and heat flows through it. 

'^ ' Assuming that fluctuation of the dipole moment, we develop an analytically solvable model for the 



problem. The effect of disorder is introduced by randomizing the positions of the dipole moments. 

We show that the disorder leads to Anderson like transition from conducting to a thermal insulating 

^ ■ state of the chain. It is shown that considered chain supports both ballistic and diffusive heat 



transports depending on the strength of the disorder. We demonstrate that nonequilibrium leads 
to the emergence of the long-range order between dipoles along the chain and make the conjecture 



^^ ' that the interplay between nonequilibrium and next to the nearest neighbour interactions results 



in the emergence of long-range correlations in low dimensional classical systems. 



1. INTRODUCTION 



The understanding of out-of-equilibrium low-dimensional systems has been a challenging 
problem over last decades. This topic covers a large variety of important problems of modern 



physics: what are the necessary conditions for the observation of the Fourier law 
to achieve and manipulate directed transport in systems with Brownian motion p, 



gain a useful work in nonequilibrium 



anspo 

m 



u 



4]; how 



6[ ; how to 



how can we control the energy transport in one- 



and two-dimensional assemblies of large organic molecules with high dipole moment arranged 



on a surface |8l-ll0[|: what is the mechanism of the transition to chaos in nonlinear chains ll|; 
what are the necessary conditions for the occurrence and existence of the temporally periodic 
and spatially localized excitations in nonlinear chains 12J ; does the unique steady state exist 
in nonlinear chains 13|. Nonequilibrium processes in low- dimensional systems are also of 
practical and techological interest because of the recent advances in nanofabrication. 

In this paper, we consider an out-of-equilibrium one- dimensional chain of particles in- 
teracting with each other via the classical dipolar potential. Our interest is twofold. First, 
we consider a heat conduction in such chains. Second, we study the emergence of new 
correlations in the system caused by a heat flow. 

We approximate a particle by a point dipole placed at a fixed position on a line. This 
approximation is valid for different real physical systems. For instance, molecules of arti- 
ficial molecular rotors contain one or several chemical groups with substantial dipole mo- 
ment. While the rest of the molecule is kept fixed on a surface these groups rotate |8|, Il0| . 
The dipole chain model is also intimately connected with ferrofluids or solid-state mag- 
netic dipoles 14J-ll6|. The point dipole approximation is valid for a the single- file water 



chain 13, UM ■ The water in narrow single wall carbon nanotubes forms a strongly ordered 

nn 

one-dimensional chain [17|, 119[. Each water molecule is connected by two hydrogen bonds 



to neighbour molecules, instead of four bonds in bulk, and also interacts with the carbon 
atoms of the nanotube. The latter interaction is weak compared to dipolar one, but ow- 
ing to the high density of the carbon atoms is not negligible and leads to the additional 



stabilization of the water in nanotube 



20 



21| . The hydrogen bonds in ID water are ener- 



getically stronger and possesses longer life-time than ones in a bulk p/7[. Both these facts 



lead to t 



pores [17|, 



le formation of the stable water molecules chain in narrow carbon nanotubes or 
22| . Moreover, equations of motion of dipole chain can be treated as a particular 




Figure 1: The sketch of the ID dipole chain. Dipoles are represented by vectors with appropriate 
magnitude and orientation. In polar coordinates the state of the dipole is described by its 
magnitude and polar angle tp. There are two heat baths that support a heat flow along a chain. 
The left bath is the "cold" one, i.e. Tl < T/j. 



case of Kuramoto model 23|, |2J]. In its turn, this model describes in a broad sense the 
synchronization in dynamical system of coupled oscillators. It is clear that this formulation 
is related to a number of phenomena from synchronization of cells in human organism to 
phase transitions or Brownian motors. The systems with dipole-dipole interaction occur not 
only in a classical physics, but also play an important role in quantum world. For instance, 
in a rapidly growing area of cold dipolar atomic gases 25|, l26|]. For example, in the high 
density limit single-file quantum dipoles can also form a lattice with a strong localization of 
atoms near lattice sites [26 1. 



The remainder of the paper is organized as follows. In Sec. 2, we describe the Langevin 
nonequilibrium dynamics for the one-dimensional chain of dipoles, theoretical approach to 
treat the disordered and the method to compute nonequilibrium correlation functions. Sec- 
tion 3 presents the results of the numerical calculations. Conclusions are given in Sec. 4. 
Some technical aspects are relegated to appendices. 



2. MATHEMATICAL FORMALISM 



A. Nonequilibrium dynamics 



Let us consider a ID chain of the A^ point dipoles with magnitude /i. The dipoles are 
rendered on a single hne with a distance a between them and are supposed to rotate in 
plane (x, y) around the axes perpendicular to the plane of the dipole. The position of the 
arbitrary dipole is characterized by a single angle </?, Figure [H Thus, our model is purely 
one-dimensional. The classical Hamiltonian of the dipole chain has the form: 



/ v^ -,9 u^ v^ /sin 09 j sin 09,- — 2 COS 09 j COS 09,- \ 



where / is the moment of inertia of the dipole, (^i is the angle between the vector of dipole 
moment of the ith dipole and the axis x, r^ is the dimensionless position (in units of a) of 
the ith dipole moment and dot stands for time derivative. 

The left and right dipoles of the chain are coupled by some mechanism to two macroscopic 



heat baths. T 
dynamics [27|, 



le energy exchange between baths and system is implemented by Langevin 



28| . The dimensionless (for units of measure see Appendix |A]) equations of 



motion have the following form: 



f 



Vi = Pi, 

N 



1 v^ 3sin(o9j + 09,-) + sin(o9j — 09,-) <- / / n .x <- / / x .x 

(2.2) 

where ru^if)-, i]R{t)) ^^^ the random forces of the left and right thermostats correspondingly. 



We assume rji^jiit) to be the Gaussian white noise 



29|. 



Left and right viscosities 7^, 7ij are 



related to noises i]L{t), i]ji{t) by the standard fluctuation-dissipation theorem: 

{V{L,R}it)mL,R}it')) = ^l{L,R}TiL,R}S{t - t'). (2.3) 

Let us assume that fluctuation of the dipoles near equilibrium positions are small enough 
to validate the expansion of sines in (12. 2 p up to the term of the first order in Lpi. This 
condition is fulfilled in the case of strong dipolar interaction and moderate temperatures of 



the thermostats. In this approximation the equations of motion become: 

Af 2 ^ m- 
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(2.4) 






E 



V'i 



tn - rj 



+ ^/i?!^) - IR^N- 



Applying the Fourier transform to both sides of (12 ■4p . we arrive at the following linear 
system of equations: 

Mc^ = T7, (2.5) 

where (p = f (^i, . . . , ^nj , t] = ir]L{uj),0, . . . ,0, r]R{^)] , 0i = exp{iut)ipi{t) dt, and 
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(2.6) 
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Now we are ready to calculate the steady-state energy current in system. To do this, let us 
consider the rate of energy change for the first dipole (connected to the left thermostat): 

d I ^ 1 ^ 

^ = 2^1 E ^1^- - 2 5^ ^^'^^ ~ ^1^^)^^ + ^^^^^^)' (2-^) 

i=2 3=2 

where ei is the energy density of the first dipole and Fij is the force that jth dipole exerts 
on first dipole. The last two terms in (12. 7p are due to the thermostat. Thus, the energy 
current that flows into the system from the left thermostat is given by: 



]{t) = -plit)-fL + pmit) = -MtfiL + Mt)vL{t). 



(2.8) 



The average current is: 



{j{t)) = {-MtflL + Mt)VL{t)). 



(2.9) 



The function v?j(t) is found from the equation (12. 5p by inverting the Fourier transform: 

V5,(u;) = {M~^)a{u)7]L{uj) + {M-'U{oo)vR{i^). (2.10) 

After some algebra we obtain the following compact expression for the steady state heat 
current in the chain: 

+cx> 
2-fLlR^T f 2wnj-l\ ( \\2 



(j) = '"^; y uj\M-^u{u)\uu, (2.11) 

— oo 

where AT = Tr — Tz,. 

The integral in (12. lip is completely defined by the poles of the (M^^)iAr element. By 
definition, poles of the M~^ are the roots of det M(a;) = equation (details of the calcu- 
lations of the integral and root finding of detM(u;) are given in Appendix [U]) . It follows 
from (12. 6p that det M(ct;) is the polynomial in a; of the order of 2A^. Therefore, according 
to fundamental theorem of algebra, there are 2A^ solutions of this equation. These solu- 
tions correspond to the frequencies of elementary excitations that carry the energy along 
the chain. Solutions of the (12.41) for free dipole chain, not connected to the thermostats, 
have form '^j{i) = e^'-* ~'^*\ where k is the wave vector and u is real. The interaction with 
thermostats leads to the occurrence of the imaginary part in u. 

The dipoles are never strictly fixed at their positions under realistic conditions. Thermal 
fluctuations and other effects are very difficult to eliminate completely. They lead to the 
fluctuations of the dipoles near their positions that brings chain into disordered state. To 
understand the role of the disorder we take into account the thermal fluctuations of dipoles 
in the direction along the line on which they are rendered. We adopt the simple scheme 
when dipoles' positions r, have Gaussian distribution 

where rjo is the position of the ith dipole in the ordered chain and a is the dispersion that 
characterizes the "strength" of the disorder. Heat current in a disordered chain is calculated 
from (12. lip by averaging over a number of realizations of the disorder. 

B. Correlation functions in nonequilibrium steady state 

The dipole orientation relaxation time is a measurable quantity in many experiments that 
allows to understand the physical nature of the processes under observation. The dipole- 



N 



dipole correlation function is defined as (/^(t) ■ /^(O)), where /^(t) = N^/j,j(t), /Xj(t) is the 

dipole moment of the ith dipole and ^{t) is total dipole moment of the chain (a polarization 
vector). It follows from the definition of the correlation function that: 



N N 

(/x(0) ■ /x(t)) = ^ ( (cos (pi{0) cos (pj{t)\ + (sin (pi{0) sin (pj{t)\) = ^ (cosAij{t) 



ij=l 



«J=1 



(2.13) 



where Ajj(t) = ipj{t) — ipi{0). It is known that if X is the Gaussian random variable than (30 1 



(exp(iX))=exp(i(X)--(X2)). 



(2.14) 



One can see that according to (12.51) (fi, i = 1,2, ... ,N, are the Gaussian random variables 
with zero mean value {ipi) = 0. Thus, we can use the rule (12.141) to calculate the average 
in dm: 



{f^{0)■^l{t))=J2^xp(-l{Al{t)) 



Z^exp ( '- + {^j{t)^i{0)) ) . 



(2.15) 



It follows from ( I2.10p that the first term in the exponent does not depend on time. The term 
{(fj{t)(fi{0)) is space cross-correlation function. It is evident from general considerations that 



r— >-+oo 



(2.16) 



Hence 



(/x(0) • ^(t)) 



i— >+oo 



E 



> > exp 



«j 



i^m) + (^f(o)) 



(2.17) 



For convenience, we subtract the rhs of this limit from (I2.15p . This will provide the dipole- 
dipole correlation function to vanish at t — > +oo. Doing this, we get: 



(Ai(0)-M(t)) = $^exp 



«j 



ivm) + i^m) 



exp((<^,(t)v9.(0)))-l 



(2.18) 



The averages in this formula immediately follows from fl2.10p : 
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+ 27i^nTRM-],{u)M-^{-u) 



duj. 



(2.19) 



The integrals are calculated by the same technique as the ones in (12.111) . 

Finally, we can write the general form of the correlation function in the linearized 
Langevin dynamics: 



N 



(MO)-M(t)) = $^K,, 



«j 



2N 



Eg., 



exp I 2_^Uij{u}k)e 
fc=i 



'ii^fci 



Kij = exp 



ivm) + ivim 



+ 1rTr 



CiN{^k)CjN{—^k) 

det'M(u;fc)detM(-u;fc) 



(2.20) 



r (,,\- — (^ T C,i{uk)Cji{-ujk) 
c^.A^fcJ 27r V"- ''det'M(a;fe)detM(-u;fc) 

det'M(a;fc) = Wjuk - ^n), 

where Uk-, k = 1,2,..., 2N are zeros of det M(co') lying in the lower half-plane of the complex 
plane u, Cij are the cofactors of the M(ci;). Coefficients G^ are calculated according the 
method described in Appendix [Bl We do not give the explicit form of Kij for brevity, but 
they can be represented in the same form as Gij. It is seen that (I2.20p decays exponentially 
only for times satisfying the condition Imukt ^ 1. 

The correlation function (I2.20p can be considerably simplified if we expand the cosine 
in (I2.13P into the power series in Ajj up to the second-order terms. In this case we have: 



iniOMt)) ^ 5^(5+(t) - lAlit)) = N'6^{t) - \ J](A?,(t)), 



«j 



(2.21) 



(A?.(t)) = ((^^(t)) + (^^^(0)) - 2(^,(t)^,(0)), 



+00 

where / 5+(t) = 1. The term (A^^.(t)) is evaluated by a direct substitution from (12.191) . 


Thus, in this simple approximation the correlation function has exponential asymptotics. 
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Figure 2: The heat conductivity of disordered dipole chain calculated according to ()2.1ip for 
different disorder magnitudes, Tl = 300 K, Tr = 400 K, Ncut = 3. On the left panel: ♦-(7 = 0, 
ordered chain; k - a = 0.02; ■ - o" = 0.05; •- a = 0.08. Right panel presents the heat conductivity 
of disordered chain in case of strong dumping, a = 0.11. The results are obtained by averaging 
over 200 realizations of disorder 



(/x(t) • /x(0)) ~ exp(— iwfct) , where Uk have the same meaning as above. Therefore, the 
longest relaxation time corresponds to the Uk with the smallest absolute value of imaginary 
part lying in the lower half-plane u. 

3. NUMERICAL RESULTS 



A. Heat conductivity and temperature profile for disordered chain 

The quantitative measure of the heat transport in media is the heat conductivity k. In 
this article we adopt the "global" definition of the heat conductivity [31] : 



k{L) 



AT' 



AT = T^ - T, 



(3.22) 



and L is the chain length. Thus, to find k we first need to find heat current. The starting 
point of calculation of heat current is (12. lip . The calculations of heat current j is straight- 
forward for ordered dipole chain. In case of disordered chain, we average a heat current over 
200 realizations of the dipole positions. The obtained dependence of the heat conductivity 
on the chain length is presented in Figure [2l It shows the transition from ballistic trans- 
port with infinite heat conductivity (for infinite chain) to diffusive transport with finite heat 
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Figure 3: Temperature profiles for different disorder strengths and corresponding localization 
lengths according to (|3.24p . The disorders (in downward direction): cr = 0.18, a = 0.03, cr = 0. 
Profiles in disordered cases were obtained by averaging over 100 disorder realizations. Dashed 
line shows position of the mean temperature T(x) = (Tl + Tr)/2. The meaning of the 
localization length (^ is revealed below in eq. ()3.24p . 

conductivity, and eventually the high level of disorder results into thermal insulating state. 
For ordered chain, a = 0, the heat conductivity is proportional to the length of the chain, 
K, r^ L. For disordered ones k(L) law deviates from the linear dependence. The observed 
change of transport regimes is an example of a very general conductor-insulator transition 



induced by a disorder [3^, |33 |. 



The disorder also affects the temperature profile of the chain. To illustrate this we 
calculated the temperature profile for different values of a. We use the kinetic definition of 
local temperature, so the temperature Tj of the zth dipole is given by 



m- 



(3.23) 



In Figure [3] we show some temperature profiles corresponding to different disorder strength 
cr. For ordered chain the temperature profile coincides with the one in harmonic chain, 
i. e. the temperature of the internal dipoles is average of the thermostats' temperatures and 
temperatures of the leftmost and rightmost dipoles are equal to the temperatures of the 
corresponding thermostats. Disorder destroys the flat temperature profile and leads to the 



10 



formation of temperature gradient. The size of the region of the chain, where it occurs, 
depends on disorder strength. Such behavior of the temperature profile could be caused 
by the localization of the elementary excitations under the influence of disorder and in this 
respect it resembles the well known Anderson transition [34 1. 

The fundamental properties of heat conduction were considered recently in concern with 
reconstruction of the Fourier's law in quantum wires |2|, |3|. Under quite general physical 
assumptions a temperature profile was shown to have the form 

AT 

^(^) = ^^ + TT 1 77Y' (3-24) 

1 + exp(-x/0 

2n n 

where x = 1, n = 1,2, . . . ,N — 1. Dubi and Di Ventra [25 dubbed ^ the "thermal 

length" because it characterizes the length-scale of the existence of a local temperature 
gradient. We approximate temperature profiles of our model by f l3.24p . In Figure [3] we see 
that increase of the disorder strength results in decrease of the ^. According to (13.241) the 
case of the ordered chain corresponds to the infinite value of ^. From Figure [3] it is seen 
that temperature gradient in the system start to develop for values of ^ of the order of a 
system size. For ^ small enough compared to the system size one observes steep temperature 
gradient near the center of a chain while the parts of the chain close to edges are thermalized 
at the temperatures of the corresponding thermostat. 

B. Nonequilibrium correlation functions and relaxation time to the steady state 

Using expression fl2.2ip . we estimate the dipole relaxation times r for different values 
of dipole moments fi. The results presented in Table [T] demonstrate inverse dependence 
t{^). It has clear physical explanation. Namely, the time of the relaxation towards an 
nonequilibrium steady state is determined by the energy transfer inside a system and higher 
magnitude of dipole moments results into more effective energy transfer between dipoles 
and, consequently, relaxation to equilibrium state runs faster. 

Two main factors influence the relaxation time in the dipole chain, see Table [H The 
first one, is the interaction strength between dipoles resulted into decreasing of a relaxation 
times, and the second one, is the chain length resulted into increasing of a relaxation time. 

To study the space correlations in chain we calculate the following correlation function 
{iPi{oo)ipj{oo)), where v^i(oo) = lim ipi{t). Details on the calculation of space correlation 



t—^+oo 
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Table I: Dipole relaxation times {Tl = 0.1 MD units, Tr = 0.3 MD units, Ncut = 3) 



/i, MD units 


1.00 


3.00 


5.00 


10.00 


15.00 


20.00 


T, ps {N = 30) 


1.25 • 10^ 


4.79 • 10^ 


3.25 • 10^ 


2.11 • 10^ 


1.72 • 10^ 


1.53 • 10^ 


r, ps {N = 50) 


6.03 • 10^ 


2.31 • 10^ 


1.57-103 


1.01 • 10^ 


8.26 • 10^ 


7.33 • 10^ 



function are given in Appendix |Bl Note, that instead of working with a number of correlation 
functions corresponding to different combinations of indexes {i,j), we construct the function 

^ N-n 



G„ 



N-n 



(3.25) 



i=l 



that expresses a " smoothed" behavior of space correlations. 

Figure H] demonstrates emergence of new long-range correlations in nonequilibrium steady 
state, which are not present in thermal equilibrium. The difference between the equilibrium 
and nonequilibrium states reaches several orders of magnitude for dipoles lying on the dis- 
tance of five and more lattice sites away from some fixed dipole. The equilibrium correlation 
function decays fast while nonequilibrium one weakly varies on the length scale of the chain. 



One more noteworthy peculiarity of space correlations is observed. Altering the range 
of dipole-dipole interaction (by changing the cutoff distance Ncut), strongly affects space 
correlations in a system. In Figure H] we see that the magnitudes of the correlations are 
substantially different for chains with Ncut = 2^4 and Ncut > 7. In Figure H] we also show 
the difference between correlations in equilibrium and nonequilibrium states. It is seen that 
there is a long range order in nonequilibrium for Ncut 7^ 2 that disappears in equilibrium 
state. Nevertheless, there is some similarity in correlations for a nearest-neighbor interaction, 
but even in this case a presence of the heat flux makes the correlations to decay slower. For 
long chains the effect still persists. 

This long-range character of the space correlation function is due to the heat flux and 
behavior of the equilibrium correlation function supports this conjecture. The right graph 
in Figure m convincingly demonstrates that without heat flux, i.e. in an equilibrium state, 
the space correlation function decays fast with distance for the same set of parameters as 
for nonequilibrium one. The link between the space correlation function and the heat flux 
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Figure 4: The absolute values of the different space correlation functions for different values of 
the interaction radius. The chain length is A^ = 50 dipoles and fi = 30 MD units, 7 = 5.0. In left 
panel the heat flux is generated by the temperature differences of thermostats Tl = 0.1 MD units, 
Tr = 0.3 MD units: T corresponds to Ncut = 2; • corresponds to Ncut = 3; ■ corresponds to 
]^cut = 4; ♦ corresponds to Ncut = 5; • corresponds to Ncut = 7; ■ corresponds to Ncut = 9; T 
corresponds to Ncut = 20. In right panel the space correlation function in equilibrium conditions 
is presented, T = 0.2 MD units, definitions of symbols are the same. 



originates from the fact that couphng to the thermostats resuhs in the interaction between 
eigenmodes of the chain 27|. From Figure [5] one can see the differences between weak 
(7 = 0.1) and strong (7 = 5.0) couphng of the chain to thermostats. 

Finally, we consider the influence of the disorder on space correlations. Surprisingly, the 
effect of the disorder turns out to be not strictly " destructive" . In Figure |6] we see that for 
a = 0.01, 0.1 correlations are stronger than ones in an ordered chain, whereas for higher 
values of a they are weaker. This result is very intriguing because from general considerations 
disorder should break the long range order in a system. 



4. SUMMARY 

In the present paper we have conducted numerical and analytical study of the classical 
dipolar chain under out-of-equilibrium conditions. We approximated nonequlibrium dynam- 
ics of the chain by a system of linearized stochastic differential equations. All the quantities 
of interest after that were expressed solely through the elements of the one matrix M ( 12. 6p . 
We focused on two basic aspects of the chain: heat conduction and correlations. 
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Figure 5: The absolute values of space correlation function in strong and weak coupling regimes. 
The chain length is iV = 50 dipoles and /i = 30 MD units, Tl = 0.1 MD units, Tr = 0.3 MD units, 
Ncut = 9: T corresponds to 7 = 0.1; • corresponds to 7 = 5.0. 
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Figure 6: The absolute values of space correlation function for different disorder strengths. The 
chain length is A^ = 50 dipoles and fi = 30 MD units, Tl = 0.1 MD units, Tr = 0.3 MD units, 
Ncut = 9: T corresponds to u = 0.00; • corresponds to a = 0.01; ■ corresponds to a = 0.10; ♦ 
corresponds to o" = 0.15; corresponds to o" = 0.20. 



To study the heat conduction we derived a closed expression for the heat current fl2.1ip . 
It was established that ordered dipole chain supports ballistic transport and their properties 
resemble ones of the harmonic lattice, e.g. the heat current is proportional to the difference 
of the thermostats' temperatures AT and not to the temperature gradient AT/N. This 
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fact points to the violation of the Fourier law in ordered dipolar chain. Ballistic transport 
regime is destroyed by a disorder introduced by a random distribution of dipoles' positions. 
We used the simple model where positions of dipoles are imposed to be random Gaussian 
variables. The dispersion of the distribution plays a role of the disorder's " strength" . Within 
the adopted model of disorder we calculated the temperature profile and heat conductivity. 
It was observed that heat conduction undergoes the transition from ballistic to diffusive 
transport. In the diffusive regime the heat conductivity decreases as we increase the chain 
length. Diffusive regime is also characterized by the establishing of the temperature gradient 
in chain. Deformation of the temperature profile in disordered chain is in good agreement 
with the recent results on the derivation of the Fourie's law in quantum wires 2|. Similarly 
to the quantum case p] , there are two different length scales in a problem of heat conduction. 
The first one corresponds to localization length. The second one corresponds to the thermal 
length. These two lengths can be very different. In the model considered in the present 
article estimation of the localization length is complicated by the presence of the strong 
correlations between eigenmodes. We left this question for further consideration. 

We constructed the exact formula for "dipole-dipole" correlation function fl2.20p . This 
allowed us to estimate the relaxation times and to show the slowing of orientation relaxation 
as the system size increases (see Table [I]). 

The situation with spatial correlations is more subtle because they are affected by different 
factors such as thermostats' temperatures difference, dipole moments and interaction range 
(cutoff radius N^ut)- The most prominent feature of nonequilibrium is the emergence of 
the long-range correlations for Ncut > 5. It is especially important because usually the 
model with only nearest-neighbor inetraction is considered in majority of works in out-of- 
equilibrium low dimensional systems. Here we clearly demonstrated that long-range behavior 
of the correlation is caused by the combination of two factors: heat flow and "long-range" 
interaction. The nonequilibrium is essential because it leads to the coupling between heat- 
carrying modes of the systems and the eigenmodes interaction is necessary for the emergence 
of long-range structure. One can generalize this conclusion by making the conjecture that 
the emergence of long-range correlations in a the one-dimensional system is possible under 
nonequilibrium conditions when next to the nearest-neighbor interactions are included. 
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Appendix A: Units of measure adopted in the article 

• The value of the dipole moment corresponds to the one of the water in the carbon 



nanotube and is taken from the molecular dynamics simulation 35|] : n = 1.9975 D; 



Length unit: a = 2.65 A, ibid; 

Energy unit: e = ^, with the above values of the dipole moment and lattice spacing 
we get e = 2.14406 ■ 10'^^ erg; 



The unit of time: r = W — — = 3.74161 • 10~^'^s, where / is the moment of inertia of 

V /^ 
a dipole. For I we take the mean value of the three principal values of this tensor of 

the water molecule 36|, / = 3.0 ■ 10^^° f^/cm^; 



Appendix B: Space correlation function in nonequilibrium steady state 

In this section we give some useful formulas that are used for calculation of the space cor 
relation function of one-dimensional lattices 37H39| . Below we closely follow the work I 



We start with rewriting the system (12. 4p in the form of the first order differential equation: 

x = Ax + F(t), (B.l) 

where the column vector x(t) = (v5i(t), . . . ,(pN{'t),Pi{'t) ■ ■ ■PNi't))'^ and 

/ Onn Inn \ „ j- . . r /t-> o\ 

^ = h Tjj = -fLOiidji + jROiNOjN (B.2) 

\ M -r / 

where Oatat is N-hj-N zero matrix, Inn is A^-by-A^ identity matrix, M is the potential energy 
matrix, F = (0Ar,?7L(t), 0, . . . ,riji{t))^ and rjLit), r]R{t) are Gaussian white noises. Solution 
of the (JB.ll) is of the form 

t 
i^it) = exp(At)x(0) + / exp (A{t - t'))F{t') dt', (B.3) 
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where x(0) is column of the initial conditions. To find the correlations in the nonequlibrium 

steady state we have to calculate the limit lim (x(t)x'(t)), where (x(t)x'(t)) is the covari- 

t— >+oo 

ance matrix. It can be done by employing the fact that all eigenvalues of the matrix A have 
negative real parts [37 



38l | that gives lim exp(At) = 0. After this remark the evaluation 

t— >-+oo 



of the limit above is done in a straightforward manner: 



lim (x(t)xt(t)) 

t— >+oo 




exp(At)Dexp(A^t)c?t, 



(B.4) 



where 



Onn Onn 

Onn a 



D = (F(t)Ft(t')) 



A,, = 2^LTL5{t - t%,5j, + 2^nTR5{t - t^^d^ 



'jN, 



(B.5) 
Here we have used the fiuctuation-dissipation theorem (12.31) and definition of F(t) given 
above. The matrices $, Z, T represent space correlation function, mean heat fiux and 
temperature profile (diagonal elements of T) correspondingly. The correlation function can 

jy method developed in [40|. It can 



be obtained by direct calculation of the integral (IB.4 
be also helpful to rewrite (IB. 40 in the matrix form [37 

AC + CA^ = -D, 



(B.6) 



with C being (x(+oo)x^(+oo)). An exact solution of this equation was found only for the 
nearest-neighbor interaction. From a general point of view, it is a well known in a control 
theory Lyapunov matrix equation and number of algorithms were developed to solve it 
numerically 4ll443|. 



Appendix C: The inverse of the matrix M and calculation of the heat current 



In the Langevin dynamics we calculate two main quantities: heat current (12.111) and the 
correlation function (I2.19p . In both cases, the formulas contain the elements of the inverse 
matrix M^^ that, by definition, is given by 

1 



M-i 



detM 
17 



C^ 



(C.l) 



where C is the matrix of cofactors [44] and T stands for transposition. Due to the M being 
symmetric matrix the matrix of cofactors C is also symmetric. 

The integrands in fl2.19p . (12. lip are the rational functions and the order of the polynomial 
in the denominator is 2N (see comments at the end of the Section 12 AD . From general 
considerations it follows that there should be no zeros of the det M(u;) on the real axis 
because if there were even one the integral in (12.111) became ambiguous while the heat 
current is physically observable quantity and must be defined unambiguously. 

Now we are ready to calculate the integral 

+ 00 

(j) = ^^^g^ / uj\M~%N{u:)?du. (C.2) 



From the residue theorem it immediately follows [45| : 



0-) = 2^.'^^^^^^ V Res (^2Cm{i^)Cm{-^) ] ^^.3) 



111 ^-^u>=^,\ detM(u;)detM(-u; 

where cjj satisfy the equation 

det M(wi) = 0, (C.4) 

in the half-plane Imw < and Res states for the residue at the point w = Wj. The 
determinant can be represented in the form of a product: 

detM(u;) = JJ(w-u;,). (C.5) 

i=X 

It is tacitly supposed that roots of this polynomial are simple. This assumption is based 
on the results of the numerical computations: for all considered chain lengths the roots are 
found to be simple and all lie in the lower half-plane of the cu-plane, Figure [71 Hence poles 
in the f IC.SP are simple and calculation of the sum can be done easily. To be assured we 
checked the multiplicity of the roots in Multroot package J46J yielded the same conclusion. 

At the end of the section we will show how to tackle the root-finding of the det M((X'). As 
it was already said, this polynomial has large coefficients and is of the order of 2A^. Hence, 
it is a cumbersome problem to find it's roots. Nevertheless, it can be greatly simplified by 
the following observation. 

Applying the Fourier transform to fIB.ip . we get the matrix equation: 



- iwx(a;) = Ax(a;) + F(w). (C.6) 
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Figure 7: Zeros of detM(a;) for different Ncut', • correspond to Ncut = 3, T correspond to 
Ncut = 20; N = 50, fi = 1.0 MD units 



Thus, a particular solution of the fIB.ip is 



+ 00 +00 

,. 1/" / ■ \ / \ , 1/" 6xp(— iu;t)_,, , , 

x(t) = -- / exp(-ia;t)x(a;) du = -—- / — — — — F(w) duj, 



(C.7) 



where I is 2N-hj-2N identity matrix. Comparison of this expression with the inverse Fourier 
transform (12. 5p shows that zeros of detM((X') coincide with zeros of det(iu;I -|- A) and the 
latter are equal to iXk, k = 1,2 . . . 2N, where A^ is the k-th eigenvalue of the A. Therefore, 

det M(u;fc) = 0, \/uJk = i\k,k = l,2,...,2N. (C.8) 

According to this observation the determinant of the lSA{u) in (]C.3p can be evaluated as 



2N 



detM(a;) = JJ(a;-iAfc). 



(C.9) 



fc=i 



The problem of finding eigenvalues of the nonsingular matrix is well known and can be 
implemented in a robust and reliable way. 

Now, when we know that all zeros of det M(a;) are simple and know the relation between 
them and eigenvalues of A, we can calculate the residue in fIC.Sp : 

Cati (iA„) Cati (— iA„) 



2N 



(j) = 2i7L7i?AT^(iA, 



n=l 



det'M(iA„)detM(-iA„)' 



(CIO) 



2N 



where det' M(iA„) = i JT(A„ — A^ 



k=l 
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Appendix D: Energy current in dipole chain 

We begin with the rate of the energy change in chain 

N , N 



23,128 



dE 



i=l i,j ^ ' ' -' ^ i,j i,j 



The energy change rate for the z-th dipole is: 

~J7 ^ nPi 2-^ -^«i ~ 9 2^ -^i«Pi • 



(D.2) 



Rewrite it in the following form: 



where 



Jl 



dej 
dt 



dt ^-^^ ■'' ' 






(D.3) 



(D.4) 



In steady state the ( — ^ ) = thus (j|") = (j°" ) for every index ?, where ( ) stands for 
the ensemble average. The j*", j°"* can be treated as the in, out - heat currents of ith dipole. 
In a straightforward manner it can be shown 27|, |28| that in a steady state m-currents are 
equal for all dipoles in chain: 



1 



32 



-P2F21 = J3 



1 ^. 
-P3 2^ Fsi 



1 



N-1 



JN 



:Pn~i 



E^' 



N-IJ- 



(D.5) 



Evidently, all these currents in steady state are equal to the energy flowed into the system 
from the "hot" thermostat per unit of time and flowed out into the "cold" thermostat. 
Equation ( ]D.5|) states that the average rate of work done by the ith dipole on the followed 
dipoles is equal to the rate of work done by the previous dipoles on the ith. Thus, one can 
say that the energy flow from the previous dipoles to the ith is equal to the energy flow from 
the ith dipole to the followed ones. 
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